#include "init.h"

int read_parameters( const char *szFileName,       /* name of the file */
		     double *Re,                /* reynolds number   */
		     double *UI,                /* velocity x-direction */
		     double *VI,                /* velocity y-direction */
		     double *PI,                /* pressure */
		     double *GX,                /* gravitation x-direction */
		     double *GY,                /* gravitation y-direction */
		     double *t_end,             /* end time */
		     double *xlength,           /* length of the domain x-dir.*/
		     double *ylength,           /* length of the domain y-dir.*/
		     double *dt,                /* time step */
		     double *dx,                /* length of a cell x-dir. */
		     double *dy,                /* length of a cell y-dir. */
		     int  *imax,                /* number of cells x-direction*/
		     int  *jmax,                /* number of cells y-direction*/
		     double *alpha,             /* uppwind differencing factor*/
		     double *omg,               /* relaxation factor */
		     double *tau,               /* safety factor for time step*/
		     int  *itermax,             /* max. number of iterations  */
		     				/* for pressure per time step */
		     double *eps,               /* accuracy bound for pressure*/
		     double *dt_value,          /* time for output */
		     int *wl,
		     int *wr,
		     int *wt,
		     int *wb,
		     double *dp,
		     char *problem,
		     int *N,
		     double *x1,
		     double *y1,
		     double *x2,
		     double *y2,
		     double *dt_path,
		     double *dt_streak,
		     double *dt_insert,
		     double *driving_velocity,
		     int *ppc
		    )
{
  READ_DOUBLE( szFileName, *xlength );
  READ_DOUBLE( szFileName, *ylength );

  READ_DOUBLE( szFileName, *Re    );
  READ_DOUBLE( szFileName, *t_end );
  READ_DOUBLE( szFileName, *dt    );

  READ_INT   ( szFileName, *imax );
  READ_INT   ( szFileName, *jmax );

  READ_DOUBLE( szFileName, *omg   );
  READ_DOUBLE( szFileName, *eps   );
  READ_DOUBLE( szFileName, *tau   );
  READ_DOUBLE( szFileName, *alpha );

  READ_INT   ( szFileName, *itermax );
  READ_DOUBLE( szFileName, *dt_value );

  READ_DOUBLE( szFileName, *UI );
  READ_DOUBLE( szFileName, *VI );
  READ_DOUBLE( szFileName, *GX );
  READ_DOUBLE( szFileName, *GY );
  READ_DOUBLE( szFileName, *PI );

  READ_INT   ( szFileName, *wl );
  READ_INT   ( szFileName, *wr );
  READ_INT   ( szFileName, *wt );
  READ_INT   ( szFileName, *wb );
  READ_DOUBLE( szFileName, *dp );

  READ_STRING( szFileName, problem);

  READ_INT   ( szFileName, *N );
  READ_DOUBLE( szFileName, *x1 );
  READ_DOUBLE( szFileName, *y1 );
  READ_DOUBLE( szFileName, *x2 );
  READ_DOUBLE( szFileName, *y2 );

  READ_DOUBLE( szFileName, *dt_path );
  READ_DOUBLE( szFileName, *dt_streak );
  READ_DOUBLE( szFileName, *dt_insert );
  READ_DOUBLE( szFileName, *driving_velocity );

  READ_INT   ( szFileName, *ppc );

  *dx = *xlength / (double)(*imax);
  *dy = *ylength / (double)(*jmax);

  return 1;
}

void init_uvp(
	      double UI,
	      double VI,
	      double PI,
	      int imax,
	      int jmax,
	      double **U,
	      double **V,
	      double **P
	      )
{
	init_matrix( U, 0, imax + 1, 0, jmax + 1, UI);
	init_matrix( V, 0, imax + 1, 0, jmax + 1, VI);
	init_matrix( P, 0, imax + 1, 0, jmax + 1, PI);
}

void init_flag(
		char *problem,
		int imax,
		int jmax,
		int **Flag
		)
{
	int i,j;
	int start_i, start_j;
	int offset;
	/* the thickness and start_i, start_j were given in the Karman spec. in worksheet 2 */
	int thickness;
	int height;

	#ifdef DEBUG
	FILE *fp;
	#endif



	/* initialze the boundary strips */
	/*
	for ( i = 0; i <= (imax + 1); i++) {
		Flag[i][0] = C_B;
		Flag[i][jmax + 1] = C_B;
	}
	for ( j = 0; j <= (jmax + 1); j++) {
		Flag[0][j] = C_B;
		Flag[imax + 1][j] = C_B;
	}



	*/

	/* Handle Karman geometry */
	if ( strncmp(problem, KARMAN, strlen(KARMAN)) == 0 ) {


		/* Fill the whole domain initially with fluid */
		for ( i = 0; i <= (imax + 1) ; i++) {
			for ( j = 0; j <= (jmax + 1); j++) {
				Flag[i][j] = C_F;
			}
		}

		thickness = 5;
		height = (jmax) / 5;
		/* Now put the obstacles */
		start_i = (2 * jmax) / 5;
		start_j = (2 * jmax) / 5;
		offset = start_i;

		for ( j = start_j ; j < (start_j + height); j++){
			for ( i = 0 ; i <  thickness; i++){
				Flag[i + offset][j] = C_B;
			}
			offset++;
		}
	}

	/* Handle Cavity geometry */
	if ( strncmp(problem, CAVITY, strlen(CAVITY)) == 0 ) {
		/* Fill the domain with fluid */
		for ( i = 0; i <= (imax + 1); i++) {
			for ( j = 0; j <= (jmax + 1); j++) {
				Flag[i][j] = C_F;
			}
		}

	}

	/* Plane geometry */
	if ( strncmp(problem, PLANE, strlen(PLANE)) == 0) {
		/* Put the pressure values */
		for ( j = 1; j <= jmax; j++) {
			Flag[0][j] = P_L;
			Flag[imax + 1][j] = P_R;
		}
		/* Fill the domain with fluid */
		for ( i = 1; i <= imax; i++) {
			for ( j = 1; j <= jmax; j++) {
				Flag[i][j] = C_F;
			}
		}

	}

	/* Step geometry */
	if ( strncmp(problem, STEP, strlen(STEP)) == 0) {

		/* Fill the domain with fluid */
		for ( i = 0; i <= (imax + 1); i++) {
			for ( j = 0; j <= (jmax + 1); j++) {
				Flag[i][j] = C_F;
			}
		}

		/* Put the pressure values */
		for ( j = 1; j <= jmax; j++) {
			Flag[0][j] = P_L;
			Flag[imax + 1][j] = P_R;
		}

		thickness = jmax / 2;
		height = jmax / 2;

		/* Now put the obstacles */
		for ( j = 0 ; j < height; j++){
			for ( i = 0 ; i <  thickness; i++){
				Flag[i][j] = C_B;
			}
		}

	}

	/* Nozzle */
	if ( strncmp(problem, NOZZLE, strlen(NOZZLE)) == 0 ) {
		/* Fill the domain with obstacles*/
		for ( i = 0; i <= (imax+1); i++) {
			for ( j = 0; j <= (jmax + 1); j++) {
				Flag[i][j] = C_B;
			}
		}


		/* Now put the fluid */
		thickness = 1;
		height = jmax - 1;
		for ( i = 0 ; i <= imax+1; i++){
			for ( j = thickness ; j <=  height; j++){
				Flag[i][j] = C_F;
			}
			thickness += 1;
			height -= 1;
			if (thickness > jmax/3)
				break;
		}
		for ( i =0 ; i <= (imax); i++ ) {
			for ( j = thickness; j <= (height); j++) {
				Flag[i][j] = C_F;
			}
		}


		/* Put the pressure values */
		for ( j = 1; j <= jmax; j++) {
			Flag[0][j] = P_L;
			Flag[imax + 1][j] = P_R;
		}
	}

	/* Reverse Flow geometry */
	if ( strncmp(problem, REVERSE_FLOW, strlen(REVERSE_FLOW)) == 0) {
		/* Put the pressure values */
		for ( j = 1; j <= jmax; j++) {
			Flag[0][j] = P_L;
			Flag[imax + 1][j] = P_R;
		}
		/* Fill the domain with fluid */
		for ( i = 1; i <= imax; i++) {
			for ( j = 1; j <= jmax; j++) {
				Flag[i][j] = C_F;
			}
		}

	}


	/* Now, we need to check the neighbours */
	for ( i = 1; i <= imax; i++) {
		for ( j = 1; j <= jmax ; j++) {

			if ( is_fluid(Flag[i - 1][j]) )
				Flag[i][j] |= F_W;
			if ( is_fluid(Flag[i + 1][j]) )
				Flag[i][j] |= F_E;
			if ( is_fluid(Flag[i][j - 1]) )
				Flag[i][j] |= F_S;
			if ( is_fluid(Flag[i][j + 1]) )
				Flag[i][j] |= F_N;

		}
	}

	#ifdef DEBUG
	fp = fopen("m.txt","w");
	fprintf(fp,"\n");

	for ( j = 0; j < (jmax + 1); j++) {
		for ( i = 0; i < (imax + 1); i++) {
			fprintf(fp,"%d ",Flag[i][j]);
		}
		fprintf(fp,"\n");
	}

	fclose(fp);
	#endif

}
